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ABSTRACT 

A fully explicit two-dimensional flow solver, based on a four-stage Runge-Kutta scheme, has been 
developed and utilized to predict two-dimensional viscous flow through turbomachinery cascades for 
which experimental data is available. The formulation is applied to the density averaged Navier-Stokes 
equations. Several features of the technique improve the ability of the code to predict high Reynolds 
number flows on highly stretched grids. These include a low Reynolds number compressible form of 
the k-e turbulence model, anisotropic scaling of artificial dissipation terms and locally varying timestep 
evaluation based on hyperbolic and parabolic stability considerations. Comparisons between 
computation and experiment are presented for both a supersonic and a low-subsonic compressor 
cascade. These results indicate that the code is capable of predicting steady two-dimensional viscous 
cascade flows over a wide range of Mach numbers in reasonable computation times. 


NOMENCLATURE 


Symbols Used 

a 

c 

c p 

Cf 

Cp 


speed of sound 
chord length 
specific heat at constant pressure 
skin friction coefficient 
pressure coefficient 
internal energy per unit mass 



e o 

E, F 
U,V 

J 

k 

L 

n, s 

P 

Pr 

Qi 

Q 

R 

S 

t 

T 

T 

A u 


total energy per unit mass (= e + V 2 /2) 

flux vectors 
contravariant velocity components 
Jacobian of curvilinear transformation 
turbulent kinetic energy 
turbulence length scale 
blade normal and tangential coordinates 

static pressure 
Prandtl number 

cartesian components of heat transfer rate vector 
primary transport variable vector 
residual vector 
source term vector 
pitch length 
static temperature 
turbulence intensity 


u, v 
V 

x, y 

a 

Pi’ p2 

5 ii 

5 w 8 ^> 
5i, 8 2 

y 

At 

e 


cartesian velocity components 
magnitude of total velocity 
cartesian coordinates 
exponent in artificial dissipation scaling function 

incidence and deviation angles 
Kronecker delta 

5.,^.^ central differencing operators 

displacement and momentum thickness 
specific heat ratio 
local timestep 

isotropic turbulent kinetic energy dissipation rate 


k 2> k 4 artificial dissipation constants 

jij, ji t molecular and turbulent viscosities 


v ij pressure monitoring parameter for artificial dissipation 

rj curvilinear coordinates 

p density 


°2> °4 artificial dissipation weighting functions 

cartesian components of stress tensor 
co _ total pressure loss coefficient 
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Superscripts and Subscripts Used 


c 

i>j 

1 

m 

o 

t 

v 

w 

0 


convective 

grid indices in streamwise and pitchwise directions respectively 

laminar 

cascade mean value (average of inlet and outlet quantities) 

stagnation 

turbulent 

viscous 

wall 

cascade pitchwise direction 


inlet ffeestream 
quantity scaled by metric Jacobian 
fluctuating quantity in time averaging 
fluctuating quantity in density averaging 
density averaged quantity 
time averaged quantity 


INTRODUCTION 

Computation of viscous flows by numerically solving the Navier-Stokes equations has become 
increasingly feasible due in most part to the ever increasing speed and memory of digital computers. 
State of the art CFD codes available today are capable of calculating steady 3-D viscous flows about 
entire vehicles, and even unsteady viscous flows in 3-D turbomachinery stages. However, despite the 
rapid advance towards exploiting the power of computers now available, some serious limitations of 
these codes have yet to be adequately resolved. Surely the most profound of these limitations is the lack 
of accurate, general turbulence models. Secondary to this, but of much concern, is the role of artificial 
dissipation in Navier-Stokes calculations. 

Explicit schemes, such as the Runge-Kutta methods first applied to the solution of the Euler equations 
by Jameson, Schmidt and Turkel 1 offer several appealing characteristics in application to fluid flow 
computations. Such schemes are easily vectorizable, amenable to convergence acceleration techniques, 
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and can be extended to unsteady flow computations in a straightforward manner. However, because of 
the stiffness associated with the explicit treatment of transport equations which contain large source 
terms, incorporation of higher order turbulence models, which contain such source terms, has not been 

popular in explicit flow solvers. 

Often, algebraic eddy viscosity models are used to approximate the apparent stresses in explicit codes. 
These models have little computational overhead and do not adversely affect the stability of the scheme. 
Though very useful in computing attached or slightly separated boundary layer flows, such models have 
well recognized drawbacks in the computation of complex flows where multiple length scales exist and 
where the transport of turbulent length scales is important. Though the k-e turbulence model also has 
major deficiencies [see Speziale 2 , Lakshminarayana 3 for example], it does provide the transport of 
length scale which is computed based on local fluid and turbulence properties. The model has been 
shown to provide better predictions than algebraic models for 2-D flow with adverse pressure gradient 
[Kirtley and Lakshminarayana 4 ] and for 2-D shock-boundary layer flow on curved surfaces [Degrez and 
VanDromme 5 ]. It therefore seemed worthwhile' to try to use it to provide an improved engineering 
approximation to the complex cascade flowfields investigated herein. 

Implicit flow solvers have been used for well over a decade to compute compressible turbulent flows 
using various forms of the k-e turbulence model. However, there have been only a few attempts to 
incorporate the model into an explicit solution procedure. In these cases, the stiffness problems 
associated with explicit treatment of the k-e model have been circumvented by incorporating semi- 
implicit treatment of the source terms [Liu 6 ], implementation of an algebraic inner layer model coupled to 
a high Reynolds number form of the k-e model in the outer layer [Liu 6 ], or by using wall functions to 
model, rather than to resolve, the near wall region where source terms and grid aspect ratio can be large 
[Grasso and Speziale 7 , Eliasson 8 , Holmes and Connell 9 ]. 

The use of higher order turbulence models, and the precise control of levels of artificial dissipation, can 
improve the accuracy of high Reynolds number flow computations about complex configurations. The 
main thrust of this investigation is the incorporation of a low Reynolds number compressible form of the 
k-e turbulence model into a purely explicit scheme, and the application of the technique to flows across a 
wide range of Mach numbers. In addition to this, some recently published improvements in controlling 
artificial dissipation levels in the computation of viscous flows on highly stretched grids are tested and 
incorporated. Two complex cascade flows are computed, for supersonic and low subsonic freestream 
conditions. For the supersonic cascade, isentropic blade Mach number, shock-boundary layer structure 
and wake loss profiles are compared with experimentally measured values. Pressure distribution and 
boundary layer profiles of velocity and turbulent kinetic energy are compared with data for the subsonic 


166 



cascade. The results are shown to be quite good within the accuracy of the turbulence model and 
experimental data used. 


GOVERNING EQUATIONS AND TURBULENCE MODEL 

In the present development, the density weighted averaging attributed to Favre 10 is used. This 
decomposition has advantages in flow computations with variable density [see Jones 11 ]. Specifically, 
the averaged governing equations are of simpler form and the physical interpretation of terms in the 
equations is clearer than when conventional time averaging is used. Reynolds (time) averaging, defined 
for a scalar, <j), as 


<{> = 4> + <})' , where <)) = 


lim i 
T— t 



<|> dt 

y 


( 1 ) 


is used for pressure, density, molecular stress tensor and molecular heat flux vector. Favre (density) 
averaging, defined for scalar, <J), as 


<j> = <J) + <j>" , where <J) = 

P , (2) 

is used for velocity components, internal energy, turbulent kinetic energy and turbulent energy 
dissipation rate. 


The resulting density averaged two-dimensional Navier-Stokes equations can be written in conservative 
form in generalized body fitted coordinates as : 
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The metric terms based on a standard coordinate transformation (x,y) — > (^,r|) are given by 


_an an 

ax’ ^ y a y ’ ax’ ^ y a y ’ 


J =£x% - ri x ^ y 


(5) 


The Jacobian, J, at each gridpoint is equal to l/(average area of adjacent cells). The contravariant 
velocity components are given by 


U = i; x u + % y v , V = T| x u + ri y v 


(6) 


Incorporating an eddy viscosity formulation, the effective stress tensor and the effective heat flux vector 
are given in cartesian coordinates by : 


rt tt 

*ij = “Clij - pu.u. = (pi 4- p t ; 


3ui 3uj 
3xj 3xi 
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( 7 ) 


. ax 

qi= - k a^ + pu > e 




In the derivation of Equations 4 and 7, it has been assumed that the time averaged molecular stress 
tensor and the molecular heat flux vector are equivalent to their density averaged values. This 
approximation should be a good one since, compared to turbulent diffusion, molecular diffusion is only 
significant near solid boundaries, where local Mach number and hence density fluctuations are small. 
The apparent heat flux vector has been modelled by incorporating the gradient diffusion hypothesis. The 
laminar Prandtl number, Pq, is set to 0.72 for air. The turbulent Prandtl number, Pr t , is set to the 
standard value of 0.90. 

In the present work, the density averaged k-e equations are numerically decoupled from the density 
averaged mean flow equations. Specifically, at each iteration, the four mean flow equations are updated 
using "frozen" values of eddy viscosity and turbulent kinetic energy from the previous iteration. 
Likewise, the coupled k-e equations are then updated using the "frozen" mean flow quantities just 
computed. 

Low Reynolds number forms of the compressible k-e equations can be written in the same form as 
Equation 3 where the scaled variable vectors become : 
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/•s P - p £ +33 

S=M, 

^(c 1 fiP-C 2 f 2 pe| + E 

where the production term P is given in cartesian coordinates as 



The mass averaged turbulent kinetic energy and isotropic component of turbulent kinetic energy 
dissipation rate are defined as 



8u- du i 
; VP 3xi"3xj 

P 


( 10 ) 


The eddy viscosity is obtained from 



( 11 ) 


The particular form of low Reynolds number model used in the code was originally devised by Chien 
for incompressible flow 12 . Compressible forms have been given by Coakley 13 and more recently by 
Nichols 14 . For this model, the constants and functions in Equations 8 through 1 1 are given by 

fj! = 1 - exp(--0ii5y + ) , fj = 1 , f 2 = 1 - 1 exp(-Rl/36) 



(12) 
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2|J.ik 


,£ = 


2 p .^8 


exp{-.5y + ) 


= 0.09 , Ci = 1.35 , C 2 = 1.80 , Pr k = 1.0 , Pr e = 1.3 


Following Hobson 15 , the blade normal coordinate, n, in Equations 12, is replaced by absolute distance 
from leading and trailing edge, upstream and downstream of the passage respectively. 

The transport variable, e, used in this model is the isotropic component of the dissipation rate, though 
the e equation is derived for the total dissipation rate. As discussed by Jones 16 , at high local Reynolds 
numbers, the anisotropic component of dissipation is negligible, so the model remains valid in these 
regions. Near a solid wall, however, the anisotropic dissipation component is not negligible, and the 
isotropic component, e, goes to zero. The term, t), accounts for the non-zero value of total dissipation 
near the wall, so that the model also remains valid near solid walls and retains the convenience of 
specifying the e - 0 boundary condition there. 

It should be noted, that although the k-e equations have been cast in compressible form, the modelling 
assumptions invoked here are essentially those for incompressible flow. Specifically, terms in the 
unmodelled k and e equations which contain density fluctuation terms, p', are neglected. Also, 
pressure diffusion terms are neglected. 

No thin layer approximations are made in either the mean flow or turbulence transport equations. 


NUMERICAL SOLUTION 


Discretization 

The H-grid flow solver used in the present studies incorporates a standard 4-stage Runge-Kutta scheme 
as first applied to Euler calculations by Jameson, Schmidt and Turkel 1 , 

Q 1 = Q + 1 At R(Q ) 

Q 2 =Q n + lAtR(Q 1 ) 

Q =Q n + lAtR(Q 2 ) 

Q =Q + At R(Q 3 ) . mi 
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Here, the residual, R, is defined according to 


R 


dE r 3F, 




aEv + aFvj +s 


(14) 


The scheme is fourth order accurate in time. Second order accurate central differences are used to 
discretize the derivatives in Equation 14. Viscous and source terms are evaluated prior to the first stage, 
convective terms are computed at every stage. The stability region for the scheme is shown in Figure 1. 



Figure 1. Stability region for standard 4- stage scheme. 


The curve in Figure 1 represents the contour lg(z)l = 1.0, where z is the complex fourier symbol of a 
discretized scalar convection-diffusion equation at a particular wavenumber, and g is the amplification 
factor arising from a 1-D scalar VonNeumann linear stability analysis of the given scheme applied to 
this discretized equation. 
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The differencing molecule used to discretize the flux vectors in computational space is shown in Figure 
2. Flux vectors, E and F, are computed at mid-points between nodes. For viscous fluxes, this scheme 

incorporates information from all nine points in the differencing molecule. This approach requires three 
times the storage for metric terms than if viscous fluxes were computed on the grid vertices themselves 
(13 point molecule). However, the truncation error associated with discretizing the viscous fluxes on a 
uniform cartesian mesh using a thirteen point scheme is 0(4Ax 2 ) + 0(4Ay 2 ) as compared to 0(Ax 2 ) + 
0(Ay 2 ) for the nine point molecule, consistent with the truncation error of the convective fluxes. It is 
also easier to apply periodic and wall-function boundary conditions with the more compact differencing 
molecule. 
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Figure 2. Computational molecule and index convention for the present scheme. Flux vectors 
computed at midpoints (marked x) between vertices. 

To accelerate the solution to steady state, locally varying timesteps are computed based on a linear 
stability analysis of the discretized Navier-Stokes equations. The resulting timestep specification is 
given as : 
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At = MIN [ Ate , At v ] = MIN 


r IMAG 

L |u| + |v| + aV(V^V^) + (Vri • Vr[) + |2(v ^ VrjJ 


REAL 

^ - v + {vr! ■ v-n)] + + m-.W v ^ • v n| + V(vfnHv^j] 

p \rri rr t r 3p 


( 15 ) 


Here, IMAG and REAL are input parameters corresponding to operational CFL and VonNeumann 
numbers chosen to ensure stability (see Figure 1). A similar expression is given by Martinelli 17 . The 
first term in the brackets in Equation 15 arises from the convection operators, the latter term corresponds 
to physical viscous terms. Note that for the present uncoupled approach, the turbulent kinetic energy 
does not appear in the stablility expressions. 

For large Reynolds number flows, the H-grids used must be highly stretched in the pitchwise direction, 
in order to adequately resolve near-wall gradients. Consequently, the metric terms, rj x , r|y in Equation 
15 can become very large near the J = 1 and J = NJ boundaries. Also, when using a two equation 
turbulence model, as in the present study, the eddy viscosity, |i t , can be very large near the J = 1 and 
J = NJ boundaries, upstream and downstream of the blades, where wall damping effects are negligible 
This combination of large eddy viscosity and grid metrics causes the viscous stability term to dominate 
in these regions, and it has been found that it is crucial to include the influence of these terms in 
determining a stable local timestep. For the supersonic and subsonic cascade test cases computed herein, 
Aty < At c for 40.0% and 7.5 % of the grid points at convergence, respectively. 


Even for the computation of steady, one-dimensional, inviscid flows, the use of highly stretched grids 
gives rise to significantly reduced convergence rates in explicit schemes. This affect arises due to 
characteristic propagation speed in the streamwise direction. In two dimensions an analagous situation 
arises when the computational mesh is clustered in one curvilinear coordinate direction to resolve regions 
where flowfield gradients are large. Specifically, the maximum local stable timestep in regions where 


the mesh is highly clustered in the ri direction, is inversely proportional to the metric term, V(Vt|-Vti) ) 
which can be very large. It is the nature of this inviscid effect which often allows one to use a local 
timestep based solely on inviscid considerations, for viscous flow computations. 


The k and e equations each contain non-linear production and destruction source terms which can be 
very large near solid boundaries. According to linear stability theory, such terms can also severely 
reduce convergence rates if a purely explicit scheme is used to discretize the equations. It was found that 
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by incorporating a composite viscous-inviscid timestep specification, the stability restrictions on the k-e 
solution are not much more severe than the restrictions on the mean flow equations discussed above. In 
fact, it was possible to compute high Reynolds number flows with this turbulence model using a purely 
explicit treatment in reasonable computation times. A local timestep for the k-e equations of 
approximately 1/4 of the stable mean flow timestep was satisfactory for the cascade flows computed 
herein. Converged solutions were thereby obtained in computation times approximately twice those of 
solutions using an algebraic eddy viscosity model. This is illustrated for turbulent flat plate flow 
computations presented in the Results section. 

Artificial Dissipation 

Central difference schemes applied to hyperbolic equations that do not contain any inherent dissipation 
require the addition of artificial dissipation, to damp high wave number disturbances. These 
disturbances can be introduced into linear problems through inconsistent boundary condition treatment 
or machine roundoff error. In non-linear problems, such disturbances can be introduced through 
aliasing of sub-grid scale non-linear disturbances, to lower, resolvable wave numbers. Even for viscous 
flow calculations, artificial dissipation must be introduced into the scheme because the physical viscous 
terms are only effective in damping frequencies at higher wave numbers than can be resolved on 
practical grids. 


In the present work, artificial dissipation is added to the discretized mean flow equations as 



dEc 3Fc 


+ 


U 


+ D(Q) + S 


(16) 


Here, D(Q) represents a mixed 2nd and 4th order nonconservative artificial dissipation operator similar 
to that devised by Jameson, Schmidt and Turkel 1 . 


D(Q) = D^(Q) + Dt|(Q) , 
D 5 (Q) = S 25 5 tt Q + S 45 8 Ktt Q , 
Dtj(Q) = §2Ti&Tyr|Q + S4 t|8t|T|T|T|Q . 


(17) 


The fourth order operators are included to damp high wave number errors and the second order 
operators are included to improve shock capturing. 
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As pointed out by Pulliam 18 , artificial dissipation terms should operate on physical values of the 
flowfield variables and as such must be appropriately scaled by the metric Jacobian. Artificial 
dissipation terms must also be scaled by the local timestep to ensure that the steady state solution is 
independent of the timestep. In addition to the above consistency requirements on the dissipation 
scaling, levels of dissipation should always be reduced to levels adequate to stabilize a scheme without 
altering the accuracy of the solution. For the computation of viscous flows on highly stretched grids, 
this latter requirement is a sensitive matter. 


For high Reynolds number flows, very highly stretched grids must be used to resolve body normal 
gradients in near wall regions. If the artificial dissipation terms in both the £, and rj directions are scaled 
by the local timestep, on grids which are highly stretched in the r\ direction, excessive dissipation is 
introduced in the £, direction. This effect is discussed by Caughey and Turkel 19 and Swanson and 
Turkel 20 . This excessive dissipation may reduce accuracy and convergence rates in viscous flow 
computations. A recently devised eigenvalue scaling of the artificial dissipation terms, due to 
Martinelli 17 alleviates this problem. Since the present technique is primarily used to compute viscous 
flows on highly stretched grids, anisotropic dissipation scaling factors similar to those used by Martinelli 
are incorporated. S 2 ^, S 2ll , S^, S 4t) in Equation 17 are defined 
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Atcq 
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At cfl 



(18) 


Here, At c ^, At CT1 , are timesteps corresponding to unit CFL limit for the inviscid one-dimensional 
problem in each direction, 

** |u| + aV(v£-V£) ^ |vj + aV(Vr|- Vq) _ (19) 

The choice of unit CFL scaling in the numerator of Equation 19 ensures that the steady state solution will 
be independent of the operational CFL limit used to compute local timesteps. If a = 1, Equation 18 
reduces to standard isotropic scaling. As mentioned above, this introduces excessive dissipation in the E, 
direction in regions where the grid is stretched in the r\ direction. If a = 0, the scaling becomes purely 
anisotropic. If the grid is very highly stretched in the q direction, such scaling may not provide enough 
dissipation in the \ direction, resulting in reduced convergence rates. For intermediate values of a 


176 



between 1/2 and 2/3, Martinelli 17 , Swanson and Turkel 20 and Radespiel and Swanson 21 have shown 
good convergence rates for Euler and Navier-Stokes calculations on highly clustered grids. 

Another scaling issue is important in the computation of viscous flows. All of the mean flow equations 
with the exception of the continuity equation contain physical dissipation terms. Also, near solid 
boundaries, physical dissipation terms in the energy equation are quite small, in the absence of heat 
transfer effects. However, near solid boundaries, the viscous fluxes in the momentum are quite large 
and are themselves adequate to provide smoothing. In these same regions, second and fourth 
derivatives of the transport variables can be quite large leading to large values of artificial dissipation 
there. This well recognized phenomenon [see Davis, Ni and Carter 22 and Swanson and Turkel 20 for 
instance] gives rise to very large nonphysical values of total dissipation in the near wall region. Often 
some sort of geometric decay function is used to control the levels of artificial dissipation in these 
regions to reduce the magnitude of numerical to physical smoothing to acceptable levels. In the present 
work, the scaling functions in Equation 18 are multiplied by a normalized square of the local velocity, 
V 2 /V 2 OOJ for the momentum equations. 

Following Jameson, Schmidt and Turkel 1 , the non-linear weighting functions in Equation 18 are 
determined from 


<*25 = k 2 max (v^i+1 j,v 5 ij,v 5 i.i j) 
<5^ = - MAX (0, K 4 - g 2 ^) 


( 20 ) 


where the monitoring parameters v, are normalized second derivatives of pressure, 




Pi+l,j " 2 Pi,j + Pi -1 j 


IPi+ld + 2 PiO + Pi-lJ 


( 21 ) 


and k 2 = 1/4, K4 = 1/64. Expressions similar to Equations 20 and 21 are used in the r\ direction. When 
shocks are not anticipated in the flowfield, k 2 is set equal to zero so that the artificial dissipation added 
to the mean flow is fourth order only. 


To examine the effects of the scalings given, two numerical experiments were conducted. The subsonic 
cascade described below was used as a numerical test bed. The reader is referred to the next section for 
specifics on this flow configuration. 
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The first experiment attempted to isolate the influence of eigenvalue scaling on accuracy and 
convergence. The test case was run for five thousand iterations using the following values of a in 
Equation 18 : a = 1 (standard isotropic scaling), a = 0 (purely anisotropic scaling) and a = 2/3 
(intermediate scaling). In Figure 3, the convergence histories for the three cases are plotted. 



# Iterations 


Figure 3. Convergence histories for various eigenvalue scalings. 

Scaling the dissipation anisotropically does provide an improved convergence rate for this case as 
expected. In addition, purely anisotropic scaling provides somewhat superior convergence rate than the 
weighted scaling (a = 2/3). This suggests that, for this case, such scaling does not reduce dissipation in 
the ^ direction to the point of destabilizing the solution. The influence of these scalings on accuracy was 
found to be negligible. 
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The second experiment sought to detect the influence of spurious dissipation levels in near wall 
boundaries, and to see how the proposed velocity scaling affects solution convergence and accuracy. 
When the ratio of artificial, to physical dissipation terms in Equation 16 were compared, it was found that 
for this test case, at convergence, artificial dissipation levels were as high as ten times the physical 
dissipation terms at the first several grid points adjacent to the wall! By incorporating the velocity 
scaling exactly as proposed above, it was possible to reduce the artificial to physical dissipation ratio to 
less than .01 in the near wall region, except in the immediate vicinity of the leading and trailing edges. 
The convergence rates compared very closely, but as shown in Figure 4, the converged solutions 
showed some discrepancy. 



s/c 

Figure 4. Comparison of predicted skin friction coefficient, along the suction surface of the test 
cascade, with and without velocity scaling of the artificial dissipation. 

It is clear from this figure that unnecessary levels of artificial dissipation in boundary layers can affect 
solution accuracy in practical application. 
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Both eigenvalue scaling (with a. = 2/3) and local velocity scaling were used in all computations that 
follow. 

It is worth noting, that the dissipation scaling considerations addressed here are especially important 
when a multigrid acceleration scheme is used. Careful tuning of artificial dissipation levels is crucial 
when performing explicit multigrid calculations on highly stretched grids. This is because inadequate or 
excessive dissipation can diminish the high wave number damping properties of the driving scheme 
thereby rendering multigrid acceleration less effective 23 . 

Lack of adequate grid resolution in the £, direction just upstream and downstream of the blade edges, 
causes the wall damping function, f^, in Equation 12 to be effective only over two to four grid points in 
this direction. This gives rise to very large streamwise gradients in k and £ at the leading and trailing 
edges, which in turn leads to slowly growing oscillations in these variables. It was found necessary to 
smooth these oscillations by incorporating small amounts of second order artificial dissipation in the k 
and £ equations in a manner consistent with Equations 16 and 17, with S 2 ^ = S 2ll = 0.002, S 4 ^ = S 4 ^ 

= 0 . 0 . 

Boundary and Initial Conditions 

Along blade surfaces the no-slip condition is imposed upon the velocities, pressure is extrapolated from 
adjacent grid points, and density is computed based on specified wall temperature or heat transfer rate. 
At the inlet, total pressure and total temperature are specified. For subsonic inflow, either inlet flow 
angle or pitchwise velocity are specified, and the R- characteristic is extrapolated along r\ = constant grid 
lines from the interior of the computational domain. At subsonic outflow boundaries, static pressure is 
specified and velocity components and entropy are extrapolated along rj = constant grid lines. Along 
periodic boundaries, cyclic information is used when discretizing derivatives in the tj direction. 

Constant values of k and £ are imposed at the inflow boundary based on specified ffeestream turbulence 
intensity and length scale, 



£ 
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Typically 5 the freestream length scale is set between .001 and .01 times the pitch of the blade passage. 
At the outflow boundary, values of k and e are extrapolated along p = constant grid lines from the 
interior of the computational domain. Turbulent kinetic energy and isotropic dissipation rate are set to 
zero along solid boundaries, as discussed in the turbulence model section. 

The flowfield is initialized using standard quasi- ID analysis to provide uniform initial velocity profiles 
along each E, - constant grid line. This gives rise to huge production terms in the k-e equations, and can 
cause solutions to become rapidly unstable. This problem is alleviated by running the code in laminar 
mode for a couple hundred iterations to develop a slight boundary layer, thereby reducing the size of 
these terms. 


RESULTS AND DISCUSSION 
Computational Considerations 

The code has been validated for laminar and turbulent flat plate boundary layer flows, where nearly exact 
agreement with theory and experiment were obtained. In addition, it has been applied to laminar flow 
about a circular arc bump in a channel, as well as to turbulent flow about a similar configuration with a 
heated wall using an algebraic eddy viscosity model. These two model problems had been computed by 
Chima and Johnson 24 and Davis, Ni and Carter 22 respectively. The present method yielded nearly exact 
agreement with these two sets of results. Implicit residual smoothing is available to accelerate the 
convergence to steady state of the mean flow. However, successful implementation of implicit 
smoothing for the turbululence transport equations has not been realized. For consistency, then, 
residual smoothing was not used in obtaining any of the proceeding results. 

For turbulent flow calculations, the highly vectorized code executes at 2.8 x 10 5 CPU seconds / 
(gridpoint * iteration) on the Cray Y-MP 8/32 at the Pittsburgh Supercomputer Center. When an 
algebraic eddy viscosity model is used, the execution rate is 1.7 x 10' 5 CPU seconds / (gridpoint * 
iteration). Since the same near wall resolution is needed for these two models, similar grids must be 
used. Experience with the code has shown that mean density residual converges slightly more slowly 
when using the k-e model, so the total overhead associated with using the higher order model is less 
than a factor of 2.0. 
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To illustrate the above considerations, convergence histories are presented here for the prediction of 
developing turbulent flow over a flat plate. Both the algebraic eddy viscosity model due to Baldwin and 
Lomax 25 and the present two-equation model were used. 

Both cases converged very slowly due to the extremely high aspect ratio of the grid (1.2 x 10 4 at the 
trailing edge of the plate). The convergence history for the computations is shown in Figure 5. It took 
approximately 10000 iterations for both calculations to converge to within engineering accuracy (taken to 
be a 4.5 order of magnitude drop in the RMS density residual). Note that the convergence rates are 
similar. This illustrates that it is primarily "inviscid" stablity constraints, and not the stiffness associated 
with large source terms in the turbulence transport equations, which give rise to the slower convergence 
rates which occur when explicit schemes are used to compute turbulent flows on highly stretched 
meshes. 


0.0 



# Iterations 

Figure 5. Convergence history for turbulent flat plate boundary layer calculations (solid line = Baldwin 

and Lomax, dashed line = k-e). 
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DFVLR PAV-1.5 Supersonic Compressor Cascade 


The first cascade to be investigated is the PAV-1.5 supersonic compressor cascade tested at DFVLR by 
Schreiber 26 . This pre-compression blade was designed especially to investigate shock-boundary layer 
interaction with separation. At the test freestream Mach number, a standoff leading edge shock forms, 
which gives rise to a separated shock-boundary layer interaction aft of mid chord on the suction surface 
of the adjacent passage. Though the measured absolute inlet Mach number was supersonic, the blade 
row stagger angle was high so the axial component of the inlet velocity was subsonic. This gives rise to 
the "unique incidence" condition wherein there exists a fixed relationship between inlet Mach number 
and inlet flow angle. Beyond a critical Mach number this condition exists and inlet conditions become 
independent of back pressure. This phenomena as well as the complex wave interaction field within the 
passage and shock-boundary-layer interaction provide a challenging test case for both numerical scheme 
and turbulence model. 

The computed case was experimentally tested in air at an inlet Mach number of 1 .53 and a maximum 
attainable static pressure ratio of 2.13. The measured axial velocity density ratio of 1.02 indicates that 
the flow was close to two-dimensional. The Reynolds number based on chord was 2.7 x 10 6 . The inlet 
turbulence intensity was measured using a Laser-two-focus (L2F) velocimeter to be no more than 1 %, 
which is the value used in the computations. As mentioned above, the inlet Mach number is supersonic, 
but axial velocity at the inlet to the computational domain is subsonic allowing left running characteristics 
to propagate out of the inlet plane. For this reason, subsonic inlet boundary conditions were specified : 
p D = 101325 N/m 2 , T 0 = 300 K, V 0OO = 379.5 m/s 2 . At the subsonic exit plane the backpressure, p e = 

56500 N/m 2 , was specified corresponding to the experimentally measured pressure ratio of the cascade, 
P 2 /P 1 = 2.13. 


The 129 x 100 computational mesh used was generated using Sorenson's 27 GRAPE code, modified by 
Gorski 28 to generate H-grids, and is shown in Figure 6. The blade normal grid spacing at the wall 
was prescribed as .00001 1 chord. This yielded values of y + < 1 at grid points adjacent to the walls. 
Except in the immediate vicinity of the leading and trailing edges, the suction and pressure surface 
boundary layers had at least 9 grid points with values of y + < 20. 
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Figure 6. 129 x 100 computational grid for the PAV-1.5 cascade. 

For clarity, only every other grid line is shown in both and ri directions. . 
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The convergence history for this computation is shown in Figure 7. It took approximately 6500 
iterations for this calculation to converge within engineering accuracy as measured by the invariance of 
total number of supersonic gridpoints in the field. This corresponded to approximately 39 minutes of 
CPU time on the Pittsburgh Cray. It was not possible to "cold start" the initialized flowfield at the 
specified pressure ratio, as the code became rapidly unstable when this was attempted. Rather, the back 
pressure had to be increased in a stepwise fashion with iteration (notice "jumps" in convergence at 
iteration 500, 1000, 2000), until the experimentally imposed pressure ratio could be specified at iteration 
2000. It is felt that the "unhealthy" convergence history is due in part to the highly clustered grid and 
also to the nearly choked operating condition. 
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Figure 7. Convergence history for PAV-1.5 cascade computation. 
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In Figure 8, a hand rendering of the shock wave pattern deduced from L2F measurements has been 
reproduced from Reference 28, alongside the computed shock wave pattern presented as divergence of 
velocity contours, and Mach number contours. 



1.05 



Figure 8. Shock wave pattern for PAV - 1.5 cascade, a) Divergence of velocity contours (-300 to 
-4800 by -500 [s' 1 ]), b) Mach number contours (0.9 to 1.5 by .15). In each diagram, the top two 
passages show computed contours. The bottom passage is the shock wave pattern deduced from flow 
visualization and L2F measurements, reproduced from Schreiber 26 . Stations labelled A-D correspond to 

.25, .50, .75 and .90 chord. 
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The key features of the flowfield are evident in this diagram, including the bow, lamda and passage 
shocks, in both experiment and computation, the bow shock is seen to impinge on the suction surface 
boundary layer of the adjacent passage. This gives rise to a lambda shock structure, a rapid thickening 
and separation of the boundary layer, and a Mach reflection which impinges on the pressure surface of 
the same passage. The high pressure ratio operating condition of this test case gives rise to a normal 
passage shock which impinges upstream of midchord on the pressure surface. This feature is also 
evident in both experiment and computation. The computation also shows some evidence of an oblique 
trailing edge shock, typical of supersonic compressor cascades at high operating pressure ratios. 

In Figure 9, the predicted isentropic blade surface Mach number is plotted against the experimental 
values. 



Figure 9. Isentropic blade surface Mach numbers for PAV-1.5 cascade computation. Calculated (solid 

line) and experimental values (symbols). 
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The calculation and experiment show fairly good agreement. The features labelled A, B and C in Figure 
9 correspond to local compression regions where the bow shock impinges on the suction surface, the 
Mach reflection impinges on the pressure surface and the passage shock impinges on the pressure 
surface. 

In Figure 10, the computed total pressure ratio is compared with traverse probe measurements at an axial 
location 0.09 chord downstream of the cascade exit plane. 
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Figure 10. Total pressure ratio profile 0.09 chord downstream of trailing edge for PAV-1.5 cascade 
computation. Calculated (solid line) and experimental values (symbols). 

The wake profile and loss distribution is reasonably well predicted, with the losses associated with the 
lambda shock system underpredicted. The wake centerline total pressure ratio is predicted reasonably 
well considering the difficulty in measurement at this location. It is noted that the results presented axe 
not fully grid independent. Modifications in the pitchwise grid clustering near midpassage gave rise to 
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and a 5% chord difference in the location of the passage shock. The loss distribution aft of the blade 
was hardly affected, but the blade surface Mach number distributions varied noticably. 


In Figure 11, computed velocity and turbulence intensity profiles at four locations on the suction and 
pressure surfaces are presented. 



Figure 11. Local velocity and turbulence intensity profiles at four chord locations along the suction 
(a and c) and pressure (b and d) surfaces for the PAV-1.5 cascade computation. The hash marks 
correspond to the estimated boundary layer thicknesses reported by Schreiber 26 . Refer to Figure 7 for 

the chord locations corresponding to A-D. 
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Along the suction surface, the predicted boundary layer is seen to remain quite thin for this high 
Reynolds number flow, thickening to about .01 chord at midchord. At .75 chord the boundary layer 
has separated due to the bow shock impinging at .67 chord, and the boundary layer thickness is seen to 
have rapidly increased to approximately .03 chord. At .90 chord, the separated boundary layer has 
grown to .05 chord. The turbulence intensity profile behaves in a manner consistent with a boundary 
layer separation. Namely, aft of the onset of separation the turbulent kinetic energy boundary layer 
thickens rapidly and the peak in intensity appears well away from the blade surface. Careful 
examination of Figure 1 1 (c) also shows that the turbulence intensity has been amplified well outside of 
the boundary layer. This amplification is presumably due to the influence of the shock on the normal 
stress components of the production term in the turbulent kinetic energy equation. 

The predicted boundary layers along the pressure surface are seen to remain quite thin along the entire 
blade. The influence of the passage shock at .40 chord is seen to thicken the boundary layer at .50 
chord, but the flow reattaches and the boundary layer thickness remains approximately .02 chord from 
.75 to .90 chord. Similar trends are noticed in the local turbulence intensity profiles, with the typical 
peak away from the blade just aft of separation, returning very close to the wall some distance after 
reattachment. 

The blade normal coordinate in Figure 1 1 is measured along grid lines which veer from perpendicular to 
the blade sufficiently far from the surface. The "kinks" in the velocity profiles at station B on the 
pressure and suction surfaces are caused by these curved grid lines intersecting the passage and bow 
shocks respectively. 

Predicted and measured performance parameters for this cascade, operating at the given conditions are 
presented in Table 1. 

Schreiber 26 , provided measured loss coefficients at maximum attainable cascade pressure ratio for a 
number of operating inlet Mach numbers. For comparison, the code was mn at two additional operating 
points within the envelope of the experimental tests. Figure 12 shows computed total pressure loss 
coefficients at all three operating points computed, along with the envelope of experimental loss 
coefficients. Computed values lie within the envelope of experimental values. It is noted, that Schreiber 
attributes the scatter in measured loss coefficient to variations in experimental axial velocity density ratio. 
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Figure 12. Total pressure loss coefficients at several cascade operating points, experiment 26 and 

computation (solid symbols). 


ARL-Double Circular Arc Subsonic Compressor Cascade 

The second cascade flow to be computed is the Applied Research Laboratory (ARL) double circular arc 
cascade tested at Penn State by Zierke and Deutch 29 . The computed case was tested at a negative 
incidence of 1.5 degrees. The working fluid was air at standard atmosphere with an inlet velocity of 
32.9 m/s (inlet Mach number = 0.1). The Reynolds number based on chord was 5.0 x 10 5 . Inlet 
turbulence intensity was measured at 1.8 %. The measured axial velocity ratio was measured to be 
between 0.97 and 1.03, indicating that the flow was close to two-dimensional. 
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Figure 13. 129 x 85 Computational grid for the ARL DCA cascade. 
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Grid spacing in the blade normal direction was set to .000023 chord on the blade surfaces. This yielded 
values of y + < 1 at grid points adjacent to the walls. Except in the immediate vicinity of the leading and 
trailing edges, the suction and pressure surface boundary layers had at least 1 1 grid points with values of 
y + < 20. 


It was only possible to obtain a steady solution when a coarse "preliminary" grid was used for this case. 
These coarse grid calculations overpredicted skin friction along the entire length of the suction surface, 
so the flow remain attached and a steady state solution was achieved. The more refined grid adequately 
resolved both inner layer and core flow regions yielding more accurate skin friction and boundary layer 
profiles. However, because both calculation and experiment show regions of mean flow reversal near 
the trailing edge, it was not possible to obtain a steady solution. The convergence history for this 
computation is shown in Figure 14. It took approximately 7000 iterations for this calculation to acquire 
a 4.5 order of magnitude drop in the RMS density residual. This corresponds to approximately 36 
minutes of CPU time on the Cray Y-MP. However as shown in Figure 14, the residual changes begin 
to increase and then level off. This is attributed to periodic shedding of vorticity from the aft portion of 
the suction surface. 


log 10 (residual) 
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Despite the lack of a steady state solution, the flow along the blade remained relatively unchanged after 
7000 iterations except for quasi-periodic shifts in the boundary layer velocity and turbulence intensity 
profiles. The measured flow also showed a small region of mean backflow near the trailing edge of the 
blade 29 , and for that reason was also probably somewhat unsteady. In Figure 15, comparison is made 
between computed blade surface pressure coefficient and measured values. Agreement is good along 
both blade surfaces. 



s/c 


Figure 15. Pressure coefficient for ARL DCA cascade computation. Calculated (solid line) and 

experimental values (symbols). 
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The oscillations in the pressure distributions near the leading and trailing edges in Figure 15 are caused 
by the velocity scaling of the artificial dissipation. The H-grid used gives rise to highly skewed regions 
near the relatively blunt leading and trailing edges of this configuration, causing the velocity scaling 
presented "as is" to give rise to these oscillations. Though the cascade flow is not significantly affected 
by this effect, it may be worth investigating improved scaling. 

In Figure 16 the predicted boundary layer profiles at three chordwise locations on the suction surface are 
plotted with those measured by laser doppler velocimeter. Agreement is excellent at 20 % chord and 
50% chord and reasonable at 90 % chord. 
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Figure 16. Boundary layer profiles at three chord locations along the suction surface for the ARL DCA 
cascade computation. Calculated (solid line) and experimental values (symbols). 
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Local turbulence intensity profiles are presented for three chord wise locations on the suction surface in 
Figure 17. As above, agreement between calculation and experiment is good at the first two stations, 
and reasonable in the aft portion of the blade. 



Figure 17. Local turbulence intensity profiles at three chord locations along the suction surface for the 
ARL DCA cascade computation. Calculated (solid line) and experimental values (symbols). 

Predicted and measured performance parameters for this cascade are also presented in Table 1. 
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Table 1. Comparison of Cascade Flow Parameters for Computed Cases 


Flowfield 

Parameter 


PAV-1.5 



ARLDCA 



Measured 

Computed 

Difference 

Measured 

Computed 

Difference 

Q a 

= 0.43 

0.38 

s 12% 

0.82 

0.88 

7% 

b 

CO 

0.144 

.149 

3.4% 

0.094 

0.111 

18 % 

Pi 

-1.8° 

-1.3° 

O 

d 

-1.5° 

prescribed 

NA 

P 2 

2.7° 

3.1° 

0.4° 

14.1° 

13.0° 

1.1° 

S ^ Cse Pss 

.63 

.67 

6.3 % 

.45 

none 

NA 

s / c sep DS 

= 0.40 

.40 

= 0.0 

none 

none 

NA 

AVDR 

1.02 

1 

NA 

0.97-1.03 

1 

NA 


C L = 


lift coefficient computed from 


(Lift per unit span) ± y n 
•5pmV^m 


CO 


pressure loss coefficients computed from 


Poe - Po 
•5P..V 2 .. 


CO = 


Po~ ' Po 


forARLand p °~" p ~ for PA V. 


CONCLUSIONS 

A Navier-Stokes procedure has been developed and applied to a supersonic and a low subsonic 
compressor cascade. A compressible low Reynolds number form of the k-e turbulence model was 
used. It was found in this study that : 

1) A fully explicit treatment of the turbulence transport equations is possible. The computational 
overhead associated with this treatment is reasonable. 

2) It is crucial to incorporate local timestep constraints based on stability analysis of the full viscous 
mean flow equations if the k-e model is used for an H-grid cascade configuration. 

3) The highly stretched grids needed to resolve near-wall physics warrant eigenvalue and local velocity 
scaling of artificial dissipation terms to improve accuracy and convergence rates. 


197 




4) Flowfield predictions were found to be good for a supersonic cascade and fair for a low subsonic 
cascade. 

5) Overall cascade performance parameters were well predicted for the supersonic cascade but not well 
predicted for the low subsonic cascade, due to flowfield unsteadiness and turbulence model 
shortcomings. 

Currently, several improvements and extensions to the technique are under way, including incorporation 
of multigridding, turbulence model corrections to account for streamline curvature and pressure strain, 
point implicit treatment of source terms in the turbulence transport equations, improved vectorization of 
the code and extension to three dimensions. 
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